%Generate the inputs for estimation in bootstrap for the main specification



clear

%&&& Import the inputs generated by STATA
choice=1:13;
ns=numel(choice);
beta=0.95;

textFileNameY=sprintf('matY.csv');  %the input file
textFileNameX=sprintf('matX.csv');   %the input file
textFileNameZ=sprintf('matZ.csv');   %the input file

tempY=dlmread(textFileNameY);  %matrix: [ahaid year hsa adj_vnd prechs sysid]
tempY=sortrows(tempY,[1,2]);
tempX=dlmread(textFileNameX);  %matrix: [ahaid year hsa bdtot prechs nprofit profit teach hhi totpop65 perc_mcr perc_mcd mktcat2 mktcat3 mktcat4 sysid]
tempX=sortrows(tempX,[1,2]);
tempZ=dlmread(textFileNameZ);  %matrix: [ahaid year vid hsa vmshbd1 vmshbd2 vmshbd3 vmshbd4 vlagbd1 vlagbd2 vlagbd3 vlagbd4 newvlag newvsysh below50cat below75cat below100cat prechs sysid]
tempZ=sortrows(tempZ,[1,2,3]);



%&&& Prepare inputs for the four types: SANA, SAAD, AFNA, AFAD
textFileName_mkt=sprintf('smplmkt.csv'); %!!!NOT yet done: Need to re-number for the sample markets used!!!
thisboot=sort(dlmread(textFileName_mkt))';  %row vector
reorder_mkt=1:numel(thisboot); %row vector
baseY=[];
baseX=[];
baseZ=[];
for i=1:numel(thisboot)
    % construct new Y
    onemkt_Y=tempY(tempY(:,3)==thisboot(i),:);
    onemkt_Y(:,3)=reorder_mkt(i);
    baseY=[baseY;onemkt_Y]; 
    % construct new X
    onemkt_X=tempX(tempX(:,3)==thisboot(i),:);
    onemkt_X(:,3)=reorder_mkt(i);
    baseX=[baseX;onemkt_X];    
    % construct new Z
    onemkt_Z=tempZ(tempZ(:,4)==thisboot(i),:);
    onemkt_Z(:,4)=reorder_mkt(i);
    baseZ=[baseZ;onemkt_Z];  
end

% ALSO need to include affiliated hospitals' brothers outside the selected
% markets in the bootstrap sample; also need to reorder mkt ID!!!
sysid=unique(baseY(:,end));
sysid(sysid==0)=[];
baseY_sys=[];
baseX_sys=[];
baseZ_sys=[];
for i=1:numel(sysid)
    % construct Y for affiliated hospitals' brothers outside the selected
    % markets; ALSO need to reorder sysid!!!
    onesys_Y=tempY(tempY(:,end)==sysid(i),:);
    incmkt=onesys_Y(:,3);  %columnn vector
    exclu=max((~(incmkt(:,ones(numel(thisboot),1))-thisboot(ones(size(onesys_Y,1),1),:))),[],2);  %only include bros OUTSIDE the selected markets
    onesys_Y(exclu==1,:)=[];
    %vmktff=max((~(Tchoice(:,ones(size(mktdom,2),1))-mktdom(ones(size(choice,2),1),:))),[],2);
    if isempty(onesys_Y)==0
       onesys_Y(:,3)=onesys_Y(:,3)+1000;  %re-order market ID to avoid replication of MKT ID with the new market ID
       baseY_sys=[baseY_sys;onesys_Y];        
    end
    % construct X for affiliated hospitals' brothers outside the selected markets
    onesys_X=tempX(tempX(:,end)==sysid(i),:);
    incmkt=onesys_X(:,3);  %columnn vector
    exclu=max((~(incmkt(:,ones(numel(thisboot),1))-thisboot(ones(size(onesys_X,1),1),:))),[],2);  %only include bros OUTSIDE the selected markets
    onesys_X(exclu==1,:)=[];
    if isempty(onesys_X)==0
       onesys_X(:,3)=onesys_X(:,3)+1000;  %re-order market ID to avoid replication of MKT ID with the new market ID        
       baseX_sys=[baseX_sys;onesys_X];        
    end
    % construct Z for affiliated hospitals' brothers outside the selected markets
    onesys_Z=tempZ(tempZ(:,end)==sysid(i),:);
    incmkt=onesys_Z(:,4);  %columnn vector
    exclu=max((~(incmkt(:,ones(numel(thisboot),1))-thisboot(ones(size(onesys_Z,1),1),:))),[],2);  %only include bros OUTSIDE the selected markets
    onesys_Z(exclu==1,:)=[];
    if isempty(onesys_Z)==0
       onesys_Z(:,4)=onesys_Z(:,4)+1000;  %re-order market ID to avoid replication of MKT ID with the new market ID                
       baseZ_sys=[baseZ_sys;onesys_Z];        
    end
end

% attach outside bros
baseY=[baseY;baseY_sys];
baseX=[baseX;baseX_sys];
baseZ=[baseZ;baseZ_sys];
    
% create new ahaid (first sort baseY, baseX and baseZ to make sure they are
% in the SAME order)
baseY=sortrows(baseY,[3,1,2]);
baseX=sortrows(baseX,[3,1,2]);
baseZ=sortrows(baseZ,[4,1,2,3]);
[~,~,newid]=unique(baseY(:,[3,1]),'rows');  % Creating IDs by groups in the order of [hsa, ahaid]
corr_newid_ahaid_newhsa=unique([newid,baseY(:,[1,3])],'rows');  % Crosswork between newid, ahaid, and reorder_mkt
id_baseZ=kron(newid,ones(ns,1));

% plug in new ahaid
baseY=[newid,baseY(:,2:end)];
baseX=[newid,baseX(:,2:end)];
baseZ=[id_baseZ,baseZ(:,2:end)];


% ### UPDATED on 12/03/2020: CANNOT use the vmshbd1-4 and vlag1-4 because {bed25,
% bed50, bed75} vary by bootstrap samples (thisboot)
bed25=prctile(baseX(baseX(:,end)==0,4),25);
bed50=median(baseX(baseX(:,end)==0,4));
bed75=prctile(baseX(baseX(:,end)==0,4),75);

all_bdtype=[(baseX(:,4)<=bed25&baseX(:,4)>0), (baseX(:,4)<=bed50&baseX(:,4)>bed25), (baseX(:,4)<=bed75&baseX(:,4)>bed50), (baseX(:,4)>bed75)];
baseZ(:,5:8)=repmat(baseZ(:,5)+baseZ(:,6)+baseZ(:,7)+baseZ(:,8),1,4).*kron(all_bdtype,ones(ns,1));
baseZ(:,9:12)=repmat(baseZ(:,9)+baseZ(:,10)+baseZ(:,11)+baseZ(:,12),1,4).*kron(all_bdtype,ones(ns,1));


% ### UPDATED on 12/03/2020: Re-calculate vsys because affiliated hospitals
% in one system could vary by bootstrap samples (thisboot)
sys_in_thisboot=unique(baseY(:,end));
sys_in_thisboot(sys_in_thisboot==0)=[];
update_vsysh=zeros(size(baseY,1),ns);
for year=2006:2009
    for i=1:numel(sys_in_thisboot)
        tmpsys=baseY(baseY(:,end)==sys_in_thisboot(i)&baseY(:,2)==year,5);  
        nummmb=numel(tmpsys);
        tmpsys(tmpsys==1)=[];                
        if isempty(tmpsys)==0
           occur=histc(tmpsys,choice);
           occur=reshape(occur,1,[]); 
           onevsysh=occur/sum(occur);
           update_vsysh(baseY(:,end)==sys_in_thisboot(i)&baseY(:,2)==year,:)=onevsysh((1:size(onevsysh,1))' * ones(1,nummmb), (1:size(onevsysh,2))' );
        end    
    end
end
baseZ(:,14)=reshape(update_vsysh',[],1);


% ### UPDATED on 12/06/2020: generage "regvar6_sing.raw"-"regvar9_sing.raw"
baseY=sortrows(baseY,[1,2]);
baseX=sortrows(baseX,[1,2]);

regvar=[baseX(baseX(:,end)==0,[1,2,4,5]),baseY(baseY(:,end)==0,[4,3])];   %matrix: newahaid year bdtot prechs adj_vnd hsa
for i=6:9
    textFileName_regvar=sprintf('regvar%d_sing.raw',i);
    output_regvar=regvar(regvar(:,2)==2000+i,:);
    output_regvar=sortrows(output_regvar,1);
    dlmwrite(textFileName_regvar,output_regvar,'delimiter', ',', 'precision', 8)
end

%the matrix: ahaid year hhi mktcat totpop65
ind_mktcat=baseX(:,13:15);
mktcat=zeros(size(baseX,1),1);
for i=1:size(baseX,1)
    if sum(ind_mktcat(i,:))>0
        mktcat(i)=find(ind_mktcat(i,:)>0);        
    end    
end
textFileName_mktcat=sprintf('mktvar_pophhi.raw');
dlmwrite(textFileName_mktcat,[baseX(baseX(:,end)==0,[1,2,9]),mktcat(baseX(:,end)==0,:),baseX(baseX(:,end)==0,10)],'delimiter', ',', 'precision', 8)




